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A  description  is  given  of  a  number  of  numerical  schemes  to  solve  an 
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evolution  equation^ that "arises  when  modelling  the  propagation  of  water  waves 
in  a  channel.  The  discussion  also  includes  the  results  of  numerical  experi¬ 
ments  made  with  each  of  the  schemes.  It  is  suggested,  on  the  basis  of  these 
experiments,  that  one  of  the  schemes  may  have  (discrete)  solitary-wave 
solutions . 
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SIGNIFICANCE  AND  EXPLANATION 


There  has  been  great  interest  in  recent  years  in  certain  nonlinear, 
dispersive  systems  used  to  model  the  propagation  of  waves  in  a  variety  of 
situations.  In  the  particular  context  of  wave  propagation  in  shallow  water, 
two  model  equations  have  been  used  widely  to  describe  the  evolution  of  such 
waves  over  a  time  scale  related  to  their  amplitude.  (A  more  complete 
discussion  of  the  time  scales  involved  and  of  the  similarity  between  the 
solutions  of  the  two  equations  is  given  in  MRC  Report  #2477  by  the  same 
authors.)  Both  of  these  models  have  smooth  solitary -wave  solutions  -  isolated 
waveforms  that  travel  at  constant  speed  and  whose  shape  is  independent  of 
time.  One  of  the  models,  the  Kbrteweg-de  Vries  equation,  has  been  studied 
extensively,  both  from  a  theoretical  and  a  numerical  viewpoint.  The  present 
paper  describes  a  numerical  study  made  of  the  second  equation.  Several 
numerical  schemes,  having  second,  fourth  and  higher  order  accuracy,  are 
described  and  results  given  of  tests  made  with  a  number  of  these  schemes 
utilizing  the  exact  solutions  to  the  continuous  problem  afforded  by  the 
solitary  wave.  A  variety  of  Bubtle  tests  were  used  to  study  the  growth  of 
errors  when  using  the  discrete  approximations,  the  results  of  which  indicated 
interesting  qualitative  differences  between  the  schemes.  In  particular,  it 
would  appear  that  one,  and  only  one,  of  the  schemes  tested  had  a  (discrete) 
solitary-wave  solution.  It  would  be  interesting  to  know  if  the  system  of 
difference  equations  associated  with  this  scheme  really  admits  a  solitary-wave 
solution. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


NUMERICAL  SCHEMES  FOR  A  MODEL  FOR 
NONLINEAR  DISPERSIVE  HAVES 


1  2  3 

J.  L.  Bona  ,  W.  G.  Pritchard  and  L.  R.  Scott 


1.  Introduction 


In  this  paper  ve  examine  some  numerical  schemes  for  the  initial-value 


problem  for  the  real-valued  function  u(x ,t)  given  by 

(PI) 
(P2) 

vhere  6  >  0  and  y  >  0  are  constants,  and  g  is  a  given  function  comprising  the 
initial  datum  for  the  differential  equation*  This  problem,  which  arises  in  the 
theory  of  water  waves,  has  been  studied  in  recent  years  by  several  workers: 
Peregrine  [13]  examined  its  possible  relevance  to  the  temporal  development  of 
undular  bores;  a  mathematical  theory  for  the  problem  was  developed  by  Benjamin, 
Bona  and  Mahony  (3] ;  and,  more  recently,  the  present  authors  |Tl  have  made  a 
detailed  comparison  of  the  model  with  the  outcome  of  some  laboratory 


ut  +  ux  +  8uux  -  y~  uxxt  -  0  ,  x  €  R,  t  >  0, 


u(x,0)  »  g(x). 


St 


experiments*  The  problem  (P)  is  closely  associated  with  the  in ital- value  problem 
for  the  Kbrteweg-de  Vries  equation 

S  ♦  ux  *  Bu\  +  y'Su  XXX  ■  °*  (1-l> 
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in  that  both  equations  have  been  advocated  as  models  for  the  same  physical 
phenomena  (e.g.  see  'Benjamin  et  al .  1 3 ] ) .  Indeed,  when  the  initial  datum  g  is 
restricted  to  conform  to  that  arising  in  many  physical  applications,  it  can  be 
shown  (see  Bona,  Pritchard  and  Scott  181)  that  the  two  equations  yield  essentialy 
the  same  solution  over  a  non-trivial  time  scale.  This  latter  work  also  points 
out  some  other  qualitative  similarities  between  the  solution  of  the  two  problems 
over  longer  time  scales. 

Several  numerical  studies  of  (P),  or  closely  associated  problems,  have  been 

reported:  e.g.  see  Peregrine  ll3l ,  Wahlbin  (l6),  Eilbeck  and  McGuire  (9),  llO]  , 
Alexander  and  Morris  |2],  Bona  et  al.  (61,  l7l»  Abdul loev  et  al .  (ll  describe  the 

results  of  some  interesting  computations  for  (P),  but  no  details  are  given  of  the 
methods  employed.  Most  of  these  studies  present  the  results  of  formal  calcula¬ 
tions,  except  for  the  work  of  Wahlbin  in  which  an  analysis  is  given  of  a  Galerkin 
method  for  a  (spatially)  periodic  version  of  (P)  and  that  of  iTl  in  which  an 
analysis  is  given  of  a  finite-difference  method  for  an  initial-  and  boundary- 
value  version  of  (P).  Both  these  latter  studies  also  showed  that  a  specific 
implementation  of  the  methods  displayed  the  expected  convergence  properties  when 
the  mesh  was  refined. 

While  developing  the  numerical  method  used  in  iTl  a  number  of  finite- 
difference  schemes  for  (P)  were  also  developed  and  tested,  and  the  purpose  of  the 
present  paper  is  to  describe  some  of  the  comparisons  that  were  made  between  these 
various  schemes.  A  description  is  given  in  §2  of  the  methods  studied,  consisting 
of  a  second-order  method  and  a  number  of  fourth-order  schemes.  The  discussion 
also  indicates  how  efficient  schemes  having  arbitrary  order  accuracy  can  be 
generated.  In  §3  a  discussion  is  given  of  the  numerical  experiments,  including 
standard  convergence  studies  along  with  a  number  of  more  subtle,  subsidiary 
experiments . 


The  main  numerical  experiments  described  below  sure  related  to  special  solu¬ 
tions  to  (P)  known  as  solitary  waves*  This  one-parameter  family  of  solutions 


represents  single-crested  waves  of  elevation  and  is  given  by 
u(x,t)  a  U  sech2{oI(x  -  xQ)  -  (l  +  —  BU)t)>, 


(1.2) 


for  U  >  0  and  a  »  I~  BY^U(l  ♦  y  BU)”*-)^2,  corresponding  to  the  initial  datura 


g(x)  =  U  sech2[a(x  -  Xq)] 


(1.3) 


The  arbitrary  parameter  Xq  gives  the  location  of  the  point  of  maximum  amplitude 
of  the  solitary  wave  at  time  t  a  o.  These  solitary  waves  propagate  without 

change  of  form  at  the  steady  speed  (l  +yBU),  determined  by  their  maximum 

amplitude  U. 

An  especially  interesting  feature  to  emerge  from  our  numerical  experiments 
is  that  one  of  the  schemes  under  study  appeared  to  have  a  discrete  solitary-wave 
solution. 
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2.  Hie  Numerical  Schemes 


In  this  section  a  description  is  given  of  the  various  schemes  that  have  been 
studied.  First  we  shall  describe  several  semi-discrete,  spatial  approximations 
to  the  solution  to  (P)  and  then  go  on  to  describe  the  temporal  approximations 
that  have  been  used.  As  indicated  in  the  introduction,  Benjamin  et_  al.  t3l 
showed  that  (P)  is  a  well-posed  problem  and  deduced  regularity  properties  of  its 
solution  u  ( in  terms  of  the  regularity  of  the  corresponding  initial  datum  g) . 
They  also  derived  an  Integral  representation  for  u,  namely 


u+(x,t)  =  /  K(x,y)(u  +  0u2)(y,t)dy, 


(2.1) 


where  K(x,y)  :=  ^  y2sgn(x  -  y)e“^x“^  ,  which  representation  we  have  used  to 
generate  the  spatial  discretizations  described  here.  Hie  ideas  outlined  below 
are  closely  related  to  the  work  described  in  [t!  •  to  which  paper  repeated 
reference  will  be  made  for  some  of  the  technical  issues  that  arise. 


2.1  Spatial  discretizations 


2.1.1.  The  GEM  scheme 

The  spatial  discretizations  were  effected  first  by  truncating  the  infinite 
interval  of  integration  to  a  finite  interval  lX1,X2l  and  then  by  taking 
quadrature  approximations  of  the  integrals 


*  IP  X« 

f  t rt  ^  n..^\  / i  \  j j  f 


/  K(x,y)(u  +  ~  Su*)(y,t)dy  and  /  iTC{x,y)(u  +  0u2)(y,t)dy.  (2.2) 

X1  x 

Justification  for  the  truncation  of  the  infinite  interval  can  be  given  using 
arguments  of  the  kind  described  in  (T 1  •  Note  that  K  is  smooth  except  for  a  Jump 
discontinuity  on  the  diagonal  y  -  x  and  so,  by  splitting  the  interval  of 
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integration  at  x,  the  smoothness  of  the  integrand  on  each  of  the  subinterivals  is 
determined  entirely  by  the  smoothness  of  the  initial  datum  g  (cf.  (2.1))*  The 
quadrature  approximation  of  the  integrals  used  here  is  the  Euler-Maclaurin 
formula  truncated  at  fourth  order,  namely  the  trapezoidal  rule  on  a  uniform  mesh 
with  one  derivative  correction  at  each  of  the  end  points  of  the  ranges  of 
integration.  When  these  derivative  corrections  fall  upon  the  unknown 
(u  +  —  0u  )  they  are  approximated  by  a  centered,  second  difference.  This 
discretization  gives  a  quadrature  rule  similar  to  the  one  derived  by  Gregory  (cf. 
Goldstine  Ill])  prior  to  the  work  of  Euler  and  Madaurin,  the  only  difference 
being  that  here  derivatives  of  K  have  been  found  exactly.  A  further  simplifi¬ 
cation  can  be  made,  as  indicated  in  (7] ,  by  ignoring  certain  small  terms  arising 
at  the  extremities  and  Xg  of  the  interval  of  integration. 

These  approximations  lead  to  a  system  of  ordinary  differential  equations  for 
functions  u^(t) ,  where  u^  approximates  u  at  the  1th  quadrature  point  (i.e», 
uA(t)  -  u( iAx ,t ) ) .  Here  i  *  Slt  Hj+1,...,»2,  with  :*  Xj/Ax,  :»  Xg/Ax, 
and  Ax  denotes  the  mesh  size.  These  equations  comprise  a  semi-discrete 
approximation  to  (2.1),  taking  the  form 

u^t)  =  Fi((u  +  |  3uou)(t)),  i  *  (2.3) 

where  u  =  (uH  ,... ,u»  ) ,  the  symbol  u®u  is  defined  by 
1  ”2 

(uou^  :=  Uj2,  (2.U) 


and  F. ( v)  =  F.  (v„  ,...,v  )  is  given  by 

1  "2 


^(v)  :=  Ax{  l  K(iAx,JAx)Vj  -  Y2(vi+1  -  v^)} 


(2.5) 


with  the  understanding  that  K(x,x)  =  0  for  all  x,  and  Vj|^  *  vNg+l  =  0  whenever 


these  expressions  appear.  ‘More  complete  details  of  the  derivation  of  these 
formulae  can  be  found  in  (7 ] •  We  shall  refer  to  the  spatial  discretization 
(2.3)-(2.5)»  as  well  as  the  simplification  to  be  described  below,  as  the  Gregory- 
Euler-Maclaurin  (GEM)  scheme. 

As  it  stands  the  method  (2.3)-(2.5)  involves  a  discrete  convolution  to 
calculate  F  »  (Fj^ .... ,Fjj^)  and  is  therefore  not  a  very  efficient  procedure. 

However,  this  may  be  overcome  in  the  following  way.  Define  a  second-difference 


operator  such  that 


(Dv^  :*  v±  -  (v1+1  -  2vt  +  v^J/te^  -  2  +  e"yAx), 


(2.6) 


which  we  write  in  the  form 


avi  +  b<vi+l  +  Vi-i). 


so  that  a  *  1  -  2b  and  -b-1  -*  (2sinh(^yAx))2  (*  (yAx)2  +  O(yAx)1*).  We  want  to 
apply  D2  to  the  convolution  term  of  (2.5),  and  it  is  therefore  convenient  to 
split  F  into  two  parts,  namely 

F^t)  =  :  Fi1(v)  +  F12(t), 


1  2  12 
where  FA  is  the  convolution  term  and  F^  (v)  :«  -  y  Ax(v^+1  - 

for  <  i  <  N2>  a  straightforward  calculation  gives 

(D2F1(v))i  *|by2Ax(vi+1  -  v1_1)  =  -12bFi2(v). 


v4  ).  Then, 


(2.7) 


If  again  we  ignore  terms  involving  points  outside  the  interval  [X1,X21,  e 
approximation  to  F*  (under  suitable  assumptions  on  v)  is  given  by  the  solution 
f(v)  to  the  tridiagonal  system  of  equations 


Thus,  it  is  more  efficient  computationally  to  use  the  semidiscrete  scheme 


u(t)  *  P((u  +  ^  Buou)(t)), 
u^O)  »  g(iAx),  i  - 


(2.9) 

(GQl) 


^  v 

where  F  ■  f  +  r  9  and  for  which  F  can  he  calculated  in  0(1*2  -'%)  operations  by 
solving  the  tridiagonal  system  (2.8).  Using  the  same  methods  as  those  described 
in  It! ,  it  can  be  shown  that  (2.9)  has  an  accuracy  of 


.  -r(X.-X- ) 

0(Ax4  +  e  2  1  ), 


(2.10) 


where  r  is  a  positive  constant. 

Although  (2.9)  does  not  appear,  super ficially,  to  be  a  standard 
discretization  it  can,  nevertheless ,  be  viewed  as  a  finite-difference 
approximation  to  (2.1)  (or  Pi)  written  in  the  form 


(1  -  y~\2\  *  -  3x(u  +  0u<:)« 

n  p 

Tb  see  this  define  difference  operators  DA  and  Dqc  by 


1  „..2, 


(Pi,  bis) 


(D1"?)*  :»  (v1+1  -  vi_1)/2Ax, 


(D02v)i  :»  -(v^  -  2vi  +  v^J/Ax* 


(2.11) 


and  a  parameter  <  by 


k  :=  -bAx2  =  (Ax/2sinh{-^  yAx))2.  (2.12) 

Then,  after  multiplying  the  (GEM)  scheme  by  D2  =  I  +  xDq2,  and  using  the 
definition  (2.8)  of  f,  it  follows  that 


(I  ♦  kD02)A  -  -y2!(tc  +  Y2*x2)I  +  ^2 


where  :*  < 


ui  +  |  eui2,  hx  <  i  <  n2. 


,  otherwise , 


(r  '3) 

jtM) 
(  ♦  ) 


and  Uj  :»  (u)it  <  i  <  N2,  u^  :a  0  otherwise. 
Note  that 


x  »  y"2(l  -  ^YAx)2  +  0((yAx)U)l  ,  (2.15) 

which,  together  with  (2.13)  can  be  used  to  generate  other  discrete  methods  for 
(PI). 


-8- 


f-  • .  Vs1 

-*  v  Vi**  v  V  * 

L?  x.  ^ *  a.  * 


’.••V./'. ^ ^ .---•wv-v-v-vl 


2.1.1  A  second-order  method. 

At  second  order  the  (GEM)  scheme  agrees  with  the  second-order  centered- 
difference  approximation  to  (Pi),  namely 

(I  +  t'2Dq2)  u  =  -DXt,  (2.16)(CD) 

with  v  defined  as  in  (2.14). 

This  kind  of  spatial  differencing  has  often  teen  used  (e.g.  see  Eilbeck  and 
^Guire  [9] ) ,  though  sometimes  the  nonlinear  term  is  not  cast  in  the 
'conservative*  form  used  here. 

2.1.3  Another  fourth-order  method 

Keeping  terms  in  (2.13)  only  to  fourth  order  yields  an  approximation  to  (Pi) 
of  Stormer-Numerov  type  (cf.  Stoer  &  Bulirsch  Il5))»  namely 

(I  +  y“2(1  -  ^  (YAx)2)DQ2)i  =  -  (I  +  Ax2^2^,  (2.17) (SN) 

and  again  v  is  defined  as  in  (2.14). 

As  with  the  (GEM)  scheme,  this  method  can  be  shown  to  satisfy  the  error 
bound  (2.10),  some  specific  tests  of  which  are  described  below. 

2.1.4  Remarks  .* 

(i)  Although  the  (GEM)  scheme  ((2.9)  or  (2.13))  has  only  fourth-order 
spatial  accuracy,  schemes  of  arbitrary-order  accuracy  can  be  derived  in  a  similar 
way  by  retaining  the  required  number  of  terms  in  the  Euler-Maclaurin  formula 
(with  the  appropriate  derivative  corrections  being  replaced  by  differences).  For 
these  higher-order  schemes  the  system  of  equations  corresponding  to  (2.8)  is  not 
altered,  the  only  change  in  the  scheme  being  that  F  is  of  the  form  F  «  f  +  F^, 
where  F^  incorporates  the  higher-order  derivative  corrections;  the  orignal  F^, 
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however,  remains  on  the  right-hand  side  of  (2.8).  This  occurs  because  the 

difference  operator  (*  I  +  <Dq^)  is  an  infinite-order  approximation  to 
—2  2 

(l  -  y  3  ),  in  the  sense  that 

d2(i  -  Y“2ax2)"1ax^  =  kdS, 

for  all  sufficiently  regular  functions  Thus  the  term  f  in  the  definition  of  F 
is  infinite-order  accurate,  and  it  remains  only  to  determine  F^  to  the  desired 
accuracy. 

(ii)  It  is  more  efficient  to  comptute  the  (GEM)  scheme  in  the  form  (2.8)- 
(2.9)  rather  than  in  the  form  (2.13):  with  the  fourth-order  scheme,  for  example, 
the  latter  arrangement  requires  the  calculation  of  a  penta-diagonal  approximation 
to  3X,  whereas  the  former  involves  only  D^.  In  general,  one  could  envisage  a 
variable-order  method  where  F®  is  calculated  to  different  orders  of  accuracy  in 
different  parts  of  the  domain  (depending,  say,  on  some  local  estimation  of  the 
spatial  errors) . . 

(iii)  the  generalisation  of  the  (GEM)  scheme  to  obtain  higher-order  methods 
may  appear  somewhat  academic,  but  the  use  of  the  fourth-order  scheme  in  [7l  (in 
modelling  a  laboratory  experiment)  placed  a  considerable  burden  on  the  data 
sampling  to  ensure  the  desired  accuracy  of  the  numerical  solutions.  Similarly, 
in  another  study  cbncerning  the  interaction  properties  of  two  solitary-wave 
solutions  of  the  family  (l.l)  (see  Bona  et  al.  [6]),  the  implementation  of  a  more 
accurate  scheme  would  have  been  beneficial.  At  the  outset  of  each  of  these 
projects  the  fourth-order  scheme  seemed  to  be  more  than  adequate  but,  in 
retrospect,  we  should  have  considered  more  seriously  the  relative  efficiency  of 
the  higher-order  schemes. 

(iv)  The  above  methods  can  readily  be  adapted  to  solve  (Pi)  posed  on  some 
fixed  interval  l ,X21  >  subject  to  the  initial  condition  that  u(x,0)  =  g(x)  for 
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Xi  <  x  <  Xg  and  the  boundary  conditions  u(x^,t)  =  h^(t)  for  t  >  0,  i  =  1,2,  where 
g,h^  and  h2  are  given  functions.  (Theory  relating  to  this  initial-  and  two-point 
boundary- value  problem  has  been  provided  by  Showalter  llU]  and  Bona  and  Doug all 8 
($)•  In  case  *  -•  or  X2  ■  +*»,  the  condition  u(x^,t)  *  h^(t)  may  be  replaced 
by  a  growth  condition,  as  in  Bona  and  Bryant  141.)  These  methods  may  also  be 
used  to  handle  the  periodic  initial-value  problem  in  which  the  initial  datum  g 
and  the  solution  u  are  both  required  to  be  periodic  in  x  with  a  given  period.  We 
have  implemented  the  GEM  scheme  for  some  of  these  problems. 

2.2  Temporal  discretization 

All  the  spatial  discretizations  of  (P)  described  above  lead  to  a  system  of 
ord lnary-di f f erent ial  equations  of  the  form 

u(t)  *  J(u(t)) ,  t  >  0, 

u(0)  «  \P , 

where,  for  example,  u^°  =  g( iAx)  for  %  <  i  <  N2.  Moreover,  the  function 
remains  suitably  bounded  as  Ax  ♦  0,  so  that  the  problem  (2.18)  is  not  in  any  way 
st iff  for  small  values  of  Ax.  An  indication  of  why  this  is  so  is  given  by.  a  von 
Neumann-type  stability  analysis  of  a  linearized  version  of  the  problem  with 

9 

periodic  boundary  conditions.  Thus,  in  (2.l4),  set  v^  =  Uu^  for  all  i,  where  U 
is  a  constant,  and  consider  the  initial-value  problem  (P)  with  2n-periodic  data 
and  corresponding  periodic  boundary  conditions  in  space.  Then,  for  all  three  of 
the  above  spatial  discretizations,  the  resulting  has  eigenvalues  Ujj  of  the 
form 

=  iUyXk, 

where  0  <  k  <  2it/Ax  and,  for  each  k,  X^  is  real  with  |XjJ  <  C,  where  C  remains 


bounded  as  Ax  approaches  zero.  In  particular,  ve  may  take 


i  *  is 


,  for  CD(2.l6), 


,  for  G1M(2.13), 


(1-^2  (y^*)2)"1/2  »  for  SN(2.17). 


Note  that,  for  the  (SN)  scheme,  ve  must  have  Ax  <  lUUy”1  in  order  that  the 
multiple  of  D02  on  the  left-hand  side  of  (2.17)  be  positive.  For  the  full 
nonlinear  problem  the  precise  boundedness  conditions  satisfied  by  $  are  given  in 
l7l  •  It  follows  that  any  of  a  variety  of  methods  for  integrating  ordinary 
differential  equations  can  be  used  to  discretize  (2.18). 

For  the  second-order  scheme  (2.16)  it  is  natural  to  consider  a  second-order 
temporal  discretization  and,  since  stiffness  is  not  a  problem,  an  explicit  method 


can  be  used.  We  have  therefore  chosen  the  so-called  "leap-frog"  scheme 

un+1  =  un-1  +  At  (u11)  ,  n  >  1. 


(2.19) 

(LF) 


The  ’starting  value’  u^  was  obtained  from  a  step  using  the  Runge-Kutta  scheme 
described  below.  (Tn  fact,  for  the  numerical  experiments  to  be  reported  in  §3, 
three  steps  of  (RK)  vere  used  Initially  so  that  (LF)  was  used  only  for  n  >  3.) 
Hie  discretization  (2.16),  (2.19)  of  (P),  the  (CD-LF)  scheme,  has  an  error  bound 


of  the  form 


«  9  -r(X,  -  X.) 

0(Ax2  +  At2  e  2  1  1 


(2.20) 


where  r  is  a  positve  constant. 
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For  the  fourth-order  schemes  (GEM)  (2.8,  2.9)  and  (SN)  we  have  used  the 
following  fourth-order  Runge-Kutta  method,  for  n  >  0: 

(2.21) 

(RK) 

Once  several  steps  of  (RK)  have  been  calculated  it  could  prove  more  efficient  to 
switch  to  a  multi-step  scheme.  We  have,  therefore,  also  considered  the  following 
prediction-correction  scheme,  for  n  >  3: 


v1  -  un  +  |  At*(un), 
v2  -  un  +  |  AtjUv1), 
w3  «  un  ♦  AtJ&(w2), 

un+1  «  un  +  £  At(2  JU1)  +  2  *(w2)  +  *(un)  +  jtfw3)] . 


uI1+1  -  un  +  ^  At 1 55  3?  -  59*n_1  +  39  J!”"2  -  9  JE®”3!, 

(2.22) 

(MS) 

un+1  »  un  +  i  At  1 251  ZiZ1*1)  *  6US3?  -  26Wn-1  +  1067?"2  -  19^“3), 


where  ^  denotes  JE(o^)  and  u1,  u2,  and  u^  are  calculated  by  (2.21).  (Note 


that  this  prediction-correction  scheme  is  used  in  the  so-called  PECE  mode,  as 
described  by  Lambert  jl2].) 

All  four  possible  combinations  of  the  spatial  discretizations  (GEM)  or  (SN), 
coupled  with  the  temporal  discretizations  (RK)  or  (MS)  provide  a  fully  discrete 
approximation  to  (P)  satisfying  the  error  bound 


U  4  -r(X  -X  ) 
0(Ax  +  At4  +  e  ), 


(2.23) 


where  r  is  a  positive  constant. 
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The  (GEM)  spatial  discretization,  coupled  with  the  standard  fourth-order 
Adams-Bashforth-Moulton  prediction-correction  scheme  (e.g.  see  Lambert  [12])  was 
used  in  IT]  to  solve  (Pi)  posed  with  initial  and  boundary  conditions* 


3*  numerical  Experiments 


3.1  Preliminary  Definitions 

In  this  section  a  description  is  given  of  some  numerical  experiments  made 
with  our  implementations  of  the  above  methods.  All  the  experiments  to  be 
described  relate  to  initial  data  g(x)  given  by  (1.3).  This  function  is 
associated  vith  the  family  (1.2)  of  exact  solutions  to  the  problem  (P)  and 
therefore  provides  a  convenient  means  of  checking  the  convergence  properties  of 
the  various  schemes.  It  has,  in  addition,  enabled  us  to  make  a  number  of  other, 
more  refined,  studies  of  the  properties  of  the  numerical  solutions. 

It  is  standard  practice  to  determine  empirically  the  convergence  of  a  scheme 
to  test  both  its  theoretical  basis  and  the  correctness  of  its  implementation. 
Fairly  detailed  studies  of  this  kind  have  been  made  for  the  methods  under 
consideration,  but  w  give  here  only  a  sample  of  the  results  that  have  been 
obtained.  All  the  experiments  to  be  described  were  made  vith  B  »  1.3  and  v2  ■  6, 
identifications  vhlch  henceforth  will  be  assigned  without  further  reference. 
These  values  relate  to  the  aforementioned  physical  problem  of  water  waves 
propagating  in  a  uniform  channel.  (We  have,  however,  carried  out  similar  tests 
vith  other  values  of  B  sad  f .)  Thus,  the  quantity  a  appearing  in  (1.2)  is  given 
by  o*  l3U/U(l  ♦  ^■U)]1^2,  and  the  solitary  wave  of  amplitude  U  has  speed 
c  s-  (1  ♦  |u). 

The  truncation  of  the  infinite  interval  described  in  §2  was  usually 
effected,  at  t  *  0,  by  choosing  X^  and  X2  such  that 

g(Xi)U  »  e  ,  i  -  1,2.  (3.1) 

Here  e  was  chosen  empirically  so  that  the  truncation  had  negligible  influence  on 
the  results.  At  each  time  step  (or  after  a  certain  specified  number  of  time 


steps)  the  right-hand  boundary  v&s  moved  outwards  (i.e.,  Xg  was  increased)  so 
that  its  distance  from  the  'crest*  of  the  solitary  wave  did  not,  on  average, 
decrease  with  time.  Whenever  such  an  expansion  of  the  domain  was  made  the  vector 
u  was  extended  by  zero  on  its  undefined  components.  In  certain  experiments  it 
was  also  possible  to  move  the  left-hand  boundary  to  the  right  without 
influencing  appreciably  the  experiment  in  question.  (In  such  cases,  the  vector  u 
was  8 imply  truncated.)  The  changes  in  and  X2  were  made  automatically  by  the 
code,  as  follows.  After  each  M  time  steps  the  values  of  X^  and  X2  were  moved  a 
distance  CjAtM,  corresponding  to  speeds  and  C2.  A  typical  value  for  M  was  25 
and  Cx  and  C2  were  chosen  appropriately  for  specific  computations.  In  particular 
<  1  ♦  ^  U  <  Cg.  We  were  conservative  about  the  positioning  of  the  endpoints 
of  the  domain  but,  even  so,  in  some  of  the  experiments  with  very  fine  meshes  the 
errors  generated  by  the  numerical  scheme  were  so  small  that  the  positioning  of 
the  boundaries  (at  t  ■  0)  gave  rise  to  a  nonnegligible  additional  error. 
However,  this  additional  contribution  was  never  more  than  5 %  of  the  total  error. 
We  shall  not,  therefore,  report  the  values  of  X1(t)  and  X2(t) ,  but  give  only 
their  values  at  t  ■  0.  These  will  be  stated  implicitly  by  quoting  either  xq  (cf. 
(1.3)),  in  which  case  X^  «  0  and  X2  a  2x0,  or  by  quoting  the  value  of  e  in  (3.1). 

To  describe  the  convergence  studies  it  is  convenient  to  introduce  some 
definitions.  Let  the  solutions  to  the  discrete  problem  at  time  t  *  JAt,  where  J 
is  a  positive  integer,  be  denoted  by  n(iAx,t),  N^t)  <  i  <  N2(t) ,  and  let  the 
exact  solution  (1.2)  be  denoted  by  u(x,t).  Define  a  relative  difference  E 
between  the  two  functions  by 

N2  N2  ' 

E(t)  :=  {  l  lu(iAx.t)  -  n(iAx,t)]2/  l  (u(iAx,t)l2}1/2.  (3.2) 

i=Nx  i=N1 

Another  functional  of  interest  in  this  problem  is  the  difference  between  the 
amplitude  nm  of  the  discrete  solution  and  that  of  the  solution  of  the  continuous 


problem,  a  difference  we  shall  refer  to  as  the  height  error  H(t).  This  quantity 
is  defined  in  the  following  way.  Find  max{n(iAx,t)  :  N^(t)  <  i  <  Ng(t)};  let  I 
be  the  value  of  1  at  which  the  maximum  is  achieved*  (if  there  is  more  than  one 
such  1,  let  I  be  the  smallest  of  these*)  Then  interpolate  the  function  values 
n(i&x,t)  by  a  quart lc  polynomial  Q(x)  at  the  points  x  *  iAx,  for  1  *  I+k  with 

|k|  <  2  (i.e.,  at  the  five  points  centered  about  the  maximum  of  n)*  Now 

determine  the  maximum  value  of  Q  using  Newton's  method,  starting  the  iteration  at 
x  *  IAx.  (This  procedure  was  successful  in  all  cases.)  Denote  this  maximum  by 
Hjg*  Then,  finally,  define  the  height  error  to  be 

H(t)  U  -  nm(t) .  (3*3) 

The  value  of  x  for  which  Q  achieves  its  maximum,  say  x^  provides  the  possibility 
of  determining  a  phase  error  at  each  t  for  the  discrete  solution  and,  by  taking 
differences,  of  obtaining  an  average  speed  for  the  wave  crest  in  the  discrete 
solution.  This  speed  can  then  be  compared  with  the  speed  c  of  the  solitary-wave 

solution  to  the  continuous  problem  and  with  the  speed  (l  +  -^  n  )  of  a  solitary 

c  m 

wave  of  the  family  (l.l)  with  amplitude  t^. 

Finally,  knowing  the  values  of  and  xm  raises  the  possibility  of  yet 
another  comparison,  namely  the  difference  between  the  solution  n(iAx,t)  of  the 
discrete  problem  and' the  function 

?(x)  »  n  sech2{[ - 2 - ]1^2(x  -  x  )}.  (3.U) 

Vl+iiL)  “ 

d.  XXL 

The  function  C  corresponds  to  the  wave  of  the  family  (l.l)  of  amplitude  nm»  whose 
crest  is  located  at  Xjj.  Then,  analogously  to  the  definition  (3.2)  we  write 

n2  n2 

D(t)  {  l  In(iAx.t)  -  s(iAx,t)]2/  l  k(iAx,t))2}1/2,  (3-5) 

i-^  i=Nx 

which  quantity  we  shall  refer  to  as  the  shape  error. 


Since,  for  the  convergence  studies,  it  is  sufficient  to  use  a  fixed  ratio  of 
Ax  to  At,  one  of  our  early  experiments  was  to  determine  an  'optimal'  value  for 
this  ratio.  The  results  of  one  such  test  (see  table  l)  indicate  that  the  best 


accuracy  was  achieved  when  Ax  ■  At.  We  decided,  therefore,  to  fix  on  the  ratio 
Ax  At  (*:  A)  for  the  remainder  of  the  study. 


At\Ax 

0.32 

0.l6 

0.08 

0.04 

0.02 

0.32 

0.25 (-1) 

0.37(-l) 

0.36(-l) 

0.36(-l) 

0.36(-l) 

0.16 

0.12(-1) 

0.12(-2) 

0.l8(-2) 

0.l8(-2) 

0.l8(-2) 

0.08 

0.l4(-l) 

0.80(-3) 

0.53( — U) 

0.86(-4) 

0.89(-4) 

o.ou 

o.i4(-i) 

0.87(-3) 

0.51(-4) 

0.28(-5) 

0.48(-5) 

0.02 

O.lU(-l) 

0.87(-3) 

0.54(-l») 

0.33(-5) 

o.i4(-5) 

TABLE  1.  The  errors  E  obtained  at  t  *  19*2  using  the  (GEM-RK)  scheme  for  a 
solitary  vave  of  amplitude  U  *  1.  (These  experiments  had  x0  3  11.)  The 
numbers  in  parentheses  Indicate  the  exponent  of  10  multiplying  the 
preceding  numbers,  e.g.  0.25 (-l)  =  0.25  x  10"1. 


A  number  of  convergence  studies  have  been  made  using  the  fourth-order 
schemes  (GEM-RK),  (GEM-MS),  (SN-RK)  and  the  second-order  scheme  (CD-LK),  the 
tests  being  carried  out  with  solitary-wave  amplitudes  U  of  0.1  and  1.0.  A 
summary  of  the  results  of  these  experiments  is  given  in  table  2  where  the  errors 
E  at  time  t  *  30,  70  and  120  are  given.  (Recall  that  the  speed  of  the 
propagation  of  these  waves  is  (l  +  ^  U)  so  that  when  t  =  120  they  will  have 
travelled  distances  of  126  and  180  spatial  units,  respectively.)  The  rows 
labelled  'ratio'  in  this  table  give  the  ratio  of  the  numbers  above  and  below  the 
entry  and  indicate  the  ratio  by  which  the  error  decreased  when  A  was  halved;  for 


the  fourth-order  schemes  this  ratio  should  approach  16  as  A  decreases  to  zero  and 
should  approach  U  in  the  same  limit  for  the  second-order  scheme*  That  this 
apparently  was  the  case  for  our  implementations  can  be  seen  from  the  results  with 
U  ■  0*1*  With  the  larger  wave  amplitude,  U  »  1.0,  convergence  orders  con¬ 
siderably  in  excess  of  4  were  found  in  many  of  the  tests  made  with  the  fourth- 
order  schemes,  suggesting  that  still  further  mesh  refinement  was  needed  before 
the  asymptotic  convergence  rate  would  be  achieved*  Under  similar  conditions, 
however,  the  asymptotic  convergence  rate  was  apparently  realized  with  the  second- 
order  (CD-LF)  scheme*  The  convergence  properties  of  the  (GEM-MS)  scheme  with 
U  «  1*0  may  seem  to  be  somewhat  anomalous ,  but  they  sure  probably  a  consequence  of 
the  fifth-order  time  stepping  used  in  this  scheme:  with  the  meshes  employed  the 
dominant  component  of  the  error  presumably  arose  from  the  temporal  approximation 
and  the  anticipated  fourth-order  rates  would'  therefore  only  emerge  with  much 
finer  meshes  or  by  using  a  different  ratio  for  At/Ax.  The  run  times  on  a  CYBER 
175  for  the  tests  with  A  ■  0.02  were  approximately  400  s  (GEM-RK),  267  s  (GEM-MS) 
and  U00  s  (SN-RK). 


30.720 


72.320 


120.320 


30.720  72.320 


(b)  (GEM-MS) 


1^0.320 


0.32 

0.108(-1) 

0.256 (-1) 

O.U29(-l) 

- 

- 

- 

ratio 

15.8 

15.7 

15.7 

- 

- 

- 

0.16 

0.685 (-3) 

0.l63(-2) 

0.27*»(-2) 

0.2U9(-2) 

0.l4l(-l) 

0.400(-l) 

ratio 

15.9 

15.8 

15.9 

28.4 

31.3 

31.3 

0.08 

0.430(-4) 

0.103<-3) 

0.172(-3) 

0.876(-4) 

0.45l(-3) 

0.128(-2) 

ratio 

16.0 

16.1 

15.9 

26.5 

32.2 

32.4 

o.o4 

0.269(-5) 

0.6i*l(— 5 ) 

0.108(-4) 

0.330(-5) 

0.l40(-4) 

0.395(-4) 

ratio 

- 

- 

- 

21.4 

24.4 

28.8 

0.02 

- 

- 

- 

0.154(-6) 

0.574(-6) 

0.137(-5) 

(a)  (GEM-RK) 

0.32 

0.107(-1) 

0.255  (-D 

0.U29(-1) 

- 

- 

- 

ratio 

15.6 

15.6 

15.6 

- 

- 

- 

0.16 

0.686(-3) 

0.l63(-2) 

0.275 (-2) 

0.167 (-1) 

0.122 

0.361 

ratio 

15.9 

15.8 

3.5.9 

13.4 

15.9 

16.6 

0.08 

0.U32(-U) 

,O.103(-3) 

0.173(-3) 

0 . 125 ( -2 ) 

0.765 (-2) 

0.2l8(-l) 

ratio 

16.0 

16.0 

16.0 

25.2 

26.8 

27.3 

O.OU 

0.270 (-5) 

0.6U3(-5) 

0.108(-4) 

0.497(-4) 

0.285 (-3) 

0.798 (-3) 

ratio 

- 

- 

- 

27.9 

29.8 

30.2 

0.02 

- 

— 

- 

0.178(-5) 

0.956(-5) 

0.26k (-4) 

•>>> 
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U  -  0.1 

A 

u 

»  1.0 

a\t 

r  30.720 

72.320 

\ 

120.320 

30*720 

72.320 

120.320 

0.32 

0.295 (-3) 

0.UU7(-3) 

0  •  563  ( —3 ) 

- 

- 

- 

ratio 

15.9 

15.9 

16.0 

- 

- 

- 

0.16 

0.186(-U) 

0.28l(-U) 

0.35l(-U) 

0.UU0(-2) 

0.l88(-l) 

O.U8l(-l) 

ratio 

16.0 

16.0 

16.0 

21. U 

25.1 

26.9 

0.08 

0.1l6(-5) 

0.176(-5) 

0.219(-5) 

0.206(-3) 

0.750(-3) 

0.179(-2) 

ratio 

15*6 

16.0 

16.0 

19.U 

22.9 

25.2 

0.0U 

0.T33(-7) 

0.110(-6) 

0.137(-6) 

0.106(-U) 

0.327 (-U) 

0.709(-U) 

ratio 

- 

- 

- 

17.2 

18.8 

21.2 

0.02 

(e)  (SN-RK) 

0.6l6(-6) 

0.17U(-5) 

0.33U(-5) 

0.32 

0.393(-2) 

0.797(-2) 

0.113(-1) 

- 

• 

• 

ratio 

3.7- 

3.9 

U.O 

- 

- 

- 

0.16 

0.106(-2) 

0.203(-2) 

0.283(-2) 

O.1U3 

0.311 

O.U96 

ratio 

3.9 

U.O 

U.O 

U.3 

U.3 

U.2 

0.08 

0.27*»(-3) 

0.513(-3) 

0.71l(-3) 

0.332(-l) 

0.725(-l) 

0.118 

ratio 

U.l 

,u.o 

U.O 

U.l 

U.l 

U.l 

0.0U 

0.670(-4) 

0.129(-3) 

0.179(-3) 

0.8l8(-2) 

0.178(-1) 

0.290(-l) 

ratio 

- 

- 

- 

U.O 

U.O 

U.O 

0.02 

- 

- 

0.20U(-2) 

O.UUU(-2) 

0.72l(-2) 

(d)  (CD-LP) 


TABLE  2.  The  errors  E  obtained  when  approximating  solitary  waves  of  amplitude  0.1 
and  1.0.  (a)  (GEM-RK)  scheme;  (b)  (GEM-MS)  scheme;  (c)  (SN-RK)  scheme; 
(d)  (CD-LF)  scheme,  e  ■  0.1  x  10"^  for  all  these  tests.  •  An  entry  in  a 
row  labelled  'ratio'  is  the  ratio  of  the  numbers  above  and  below  that 
entry. 
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While  table  2  indicates  the  convergence  properties  of  the  schemes  with  A, 
the  graphs  in  figure  1  indicate  the  temporal  dependence  of  E  for  a  fixed  A. 
Figure  la  shows  the  error  E  for  the  various  approximations  computed  with  A  =  0.0U 
to  the  solitary  wave  having  U  *  0.1,  and  figure  lb  shows  the  approximations  found 
with  A  *  0.02  to  the  wave  having  U  *  1.0.  (Note  that  E  is  plotted  on  a 
logarithmic  scale.)  The  (SN-RK)  scheme  gave  tue  best  approximation  to  the 
smaller  solitary  wave  whereas  the  (GEM-RK)  scheme  gave,  by  and  large,  the 
smallest  errors  for  the  solitary  wave  with  U  =  1.0.  (Note  that  the  (GEM-RK)  and 
the  (GEM-MS)  schemes  gave  nearly  the  same  errors  for  the  computation  of  the 
smaller  solitary  wave,  but  for  the  larger  wave  the  errors  were  greatly  different, 
suggesting  that  the  choice  At  =  Ax  was  not  1  optimal1  for  the  (GEM-MS)  scheme.) 
Both  sets  of  graphs  show  an  Initial  phase  over  which  the  error  rapidly  increased 
and  after  which  there  was  a  slower  increase  in  E.  Much  of  the  slower  increase 
arises  because  the  'amplitude'  of  the  approximate  solutions  is,  in  general, 
different  from  that  of  the  exact  solution  and  therefore  the  two  waves,  having 
slightly  different  phase  speeds,  draw  apart. 

3.3  Further  Tests 

The  dependence  of  the  height  error  H  on  time  for  the  various  schemes  (when 

f 

A  =  0.02)  is  shown  in  figure  2a.  The  results  for  figure  2  were  obtained  using  a 
solitary  wave  of  amplitude  U  =  1.0  as  the  initial  datum.  (Note  that,  to 
normalize  the  errors,  the  ordinates  in  figure  2  have  been  scaled  by  An,  where  n 
is  the  order  of  accuracy  of  the  scheme  being  used.)  Thus  we  see  that  the 
'amplitude'  of  the  discretely-computed  waveforms  for  both  the  (GEM-RK)  and  the 
(SN-RK)  schemes  were  slightly  smaller  than  the  amplitude  of  the  exact  solution; 
moreover,  after  an  initial  'settling  out'  period,  the  amplitude  of  the 
approximate  solution  decreased  monotonically  with  time.  Both  the  (GEM-MS)  and 
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FIGURE  2.  Two  measures  of  differences  between  approximate  solutions  and  solitary^ 
wave  solutions,  for  calculations  using  initial  data  corresponding  to  a 


solitary  wave  of  amplitude  U  *  1.0  (and  c  »  0.1  *  10"^). 


(GEM-RK), 


n  *  4;  —  •  —  :  (GEM-MS),  n  *  4; - -  (SN-RK),  n  »  4; - (CD-LF), 

n  ■  2.  These  calculations  had  ^  ■  0.02.  (a)  H(t)/An  (see  3.3); 

(b)  D(t)/An  (see  3.5). 


the  ( CD-LF )  schemes ,  on  the  other  hand,  generated  waveforms  whose  amplitudes 
exceeded  that  of  the  exact  solution.  Probably  the  most  interesting  feature  of 
these  computations  is  that  the  (CD-LF)  scheme,  in  contrast  to  the  other  schemes, 
generated  a  waveform  vhich,  after  an  Initial  period,  had  an  amplitude  (or  height 
error)  independent  of  t.  We  shall  consider  this  point  in  more  detail  below. 

Another  obvious  feature  of  the  graphs  is  the  'kink'  near  t  =*  10.  This 
corresponds  to  the  time  at  which  the  crest  of  the  wave  passed  the  initial 
location  of  the  right-hand  boundary  of  the  domain,  the  point  where  the  Initial 
datum  had  been  truncated.  This  is  not  unexpected  since  the  theory  for  problem 
(P)  developed  in  [3l  showB  that  a  discontinuity  of  the  sort  generated  by  our 
truncation  procedure  does  not  propagate.  (The  larger  errors  associated  with  the 
second-order  scheme  presumably  dominated  the  truncation  effect  so  that  the  'blip' 
was  not  apparent  in  that  case.) 

The  same  effect  is  also  apparent  in  connection  with  the  ’shape  error'  D(t) , 
shown  in  figure  2b,  where  a  'spike'  can  be  seen  near  t  *  10  for  all  three  of  the 
higher-order  schemes.  Recall  that  the  shape  error  indicates  the  difference 
between  the  discretely-computed  waveform  and  a  solitary  wave  (l.l)  having  the 
same  amplitude  and  phase  location  as  that  of  the  discrete  solution.  Thus,  the 
results  of  figure  2b  suggest  that,  after  an  initial  'settling-out'  period  of 
about  30  time  unitrf,  the  shape  of  the  discretely-computed  wave  changed  only 
rather  slowly  with  time,  with  the  exception  of  the  (GEM-MS)  scheme  which  showed 
some  short-time  variation  in  D  superimposed  on  a  more  gradual,  long  term 
variation.  In  keeping  with  the  results  shown  in  figure  2a,  the  (CD-LF)  scheme 
generated  a  waveform  whose  shape  was  eventually  independent  of  t. 

Thus,  the  above  results  suggest  that  the  (CD-LF)  scheme  has  discrete 
'solitary-wave'  solutions,  whereas  none  of  the  other  schemes  would  appear  to  have 
this  property.  lb  provide  further  support  for  this  thesis  similar  tests  were 


made  with  a  variety  of  different  initial  wave  amplitudes  U.  This  new  set  of 
experiments  was  made  with  A  =  0.l6  with  the  calculations  progressing  for  8000 
time  steps  to  investigate  the  possibility  of  very  slow  temporal  changes  in  H  or 
D.  The  results  of  these  tests  are  summarised  in  figures  3  and  U,  the  height 
errors  being  given  in  figure  3  and  the  shape  errors  in  figure  U.  It  is  seen  from 
these  graphs  that  the  height  error  for  the  (CD-LF)  scheme  quickly  became  constant 
at  a  value  that  was  maintained  for  severed  thousand  time  steps*  By  contrast  the 
height  errors  for  both  the  (GEM-RK)  and  the  (SN-RK)  schemes  Increased  steadily 
with  time  (see  figure  3a  where  the  case  U  =  1.0  is  shown).  An  illustration  of 
how  nearly  constant  the  height  error  was  for  the  ( CD-LF)  scheme  is  as  follows. 
At  t  *  23.0U,  H  =  -0.01TU6  (rounded  to  significant  digits)  and  from  then  until 
t  *  1280.00  the  height  error  was  either  -0.017^6  or  -0.017U7,  with  the  fluctua¬ 
tion  in  the  fourth  digit  probably  arising  as  a  consequence  of  the  height- locating 
procedure  that  was  used  (cf.  §3.1) •  The  times  taken  for  the  steady  situation  to 
be  attained  were  found  to  depend  rather  strongly  on  the  amplitude  U  of  the 
Initial  datum,  a  feature  that  is  evident  in  figure  3,  but  which  appears  more 
obviously  in  the  graphs  of  the  shape  error  shown  in  figure  U.  Thus,  when  U  =  1.0 
it  took  only  about  U0  time  units  for  the  steady  waveform  to  be  realized  whereas 
with  U  *  0.1  it  took  at  least  1000  time  units  for  the  discrete  waveform  to  reach 
its  steady  value. 

Also  given  in  figure  ta  are  graphs  of  the  shape  error  D  for  both  the  (GEM- 
RK)  and  the  (SN-RK)  schemes,  for  the  case  U  =  1.0.  It  can  be  seen  from  these 
graphs  that  the  variations  of  D  with  time  were  quite  small  when  t  exceeded  about 
100.  Rote  that  this  does  not  mean  the  discrete  wave  nearly  has  permanent  form 
for,  as  we  saw  in  figure  3,  the  'amplitudes’  of  these  discrete  solutions 
decreased  steadily  with  time.  Rather,  the  interpretation  is  that,  for  t  >  100 
say,  the  'shape*  of  these  discrete  solutions  differed  from  a  solitary  wave  (l.l) 


FIGURE  3.  Th*  height  error  when  A  ■  0.16  for  *  variety  of  Initial  wave 

amplitude*  U.  (c  ■  0.1  a  10"^).  Rote  the  ordinate*  are  magnified  by  the 
factor  10s. 


(a) 


(b) - 


1.0,  (GEM-RK)i  „ 
J  a 


(S3-RK) 

U  ■  1.0,  (CD-LF), 
U  ■  0.T5,  (CD-LF) , 
U  ■  0.5,  (CD-LF), 
U  •  0.25,  (CD-LF), 
U  ■  0.1,  (CD-LF), 


U,  s 


n  ■  2,  s 
n  ■  2,  s 
n  •  2,  s 
n  ■  2,  s 
n  ■  2,  s 


-2,  H(l280)/An  »  3 **.86 


1, 

1. 

1, 

2, 

3, 


H(1280)/AB  »  -0.6T22; 
H(l280)/An  *  -0 . 3**80 ; 

h( 1280) /a"  »  -O.IU12. 


FIGURE  It.  the  shape  error  vtaen  A  ■  0.16  for  a  variety  of  Initial  wave  anpli 
tildes  U.  (c  ■  0.1  *  10“^.) 


(a)  —  •••  —  :  U  »  1.0, 

(OEM-RK), 

n  »  J», 

D( 1280 )/4n  *  0.8869; 

-  :  U  ■  1.0, 

(SR-RX), 

n  »  U, 

D(l280)/An  «  0.7119; 

—  •  —  :  U  ■  1.0, 

(CD-LF), 

n  ■  2, 

D(l280)/An  *  0.9**l6; 

-  :  U  »  0.T5, 

(CD-LF), 

n  ■  2, 

D(l280)/An  -  0.6263; 

-  •  u  »  0.5, 

(CD-LF), 

n  ■  2, 

D(l280)/An  *  0.3680. 

(b) - -  u  -  0.25, 

(CD-LF) , 

n  ■  2, 

D(  1280 )/An  »  0.1605; 

-  :  U  -  0.1 

(CD-LF) , 

n  ■  2, 

D(l280)/An  *  0.05371 

of  the  sane  amplitude  by  approximately  a  constant  proportion.  Thus,  although  the 
amplitudes  of  these  approximations  to  the  solitary  wave  (l.l)  decreased  steadily 
with  time  the  waveform  preserved  a  shape  fairly  close  to  that  of  a  member  of  the 
family  (l.l). 

Some  of  the  values  of  D/A**  for  the  (GEM-RK)  scheme,  from  which  the  graph  in 
figure  4a  was  drawn  are  given  in  table  3* 


99.84 

199.68 

299.52 

399.36 

499.20 

599.04 

701.44 

0.8777 

0.8812 

0.8825 

0.8838 

0.8865 

0.8860 

0.8855 

801.28 

901.12 

1000.96 

1100.80 

1200.64 

0.8841 

0.8850 

O.8829 

0.8835 

0.8795 

TABLE  3.  Some  values  of  D/A1*  for  the  (GEM-RK)  scheme,  U  »  1.0, 
e  *  0.1  x  10_T,  A  ■  0.16,  corresponding  to  the  graph 
shown  in  figure  4a. 
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k.  Conclusions 

A  description  has  been  given  of  a  number  of  numerical  schemes  to  solve  the 
initial-value  problem  (P).  The  methods  studied  have  either  second-  or  fourth- 
order  accuracy  in  both  the  space  and  time  variables ,  though  one  of  the  schemes 
has  a  straightforward  extension  to  any  desired  order  of  accuracy.  Rigorous  error 
estimates  can  be  obtained  for  all  these  schemes  using  the  methods  described  in 
IT  1  •  A  description  is  also  given  of  some  numerical  experiments  made  with  these 
schemes  based  on  the  solitary-wave  solution  (1.2)  of  (P).  The  experiments,  which 
included  both  a  standard  convergence  study  and  other  special  tests,  revealed  some 
subtle  differences  between  the  errors  for  the  various  schemes.  So,  for  example, 
the  multi-step  scheme  appeared  to  introduce  an  oscillatory  component  to  the 
solution,  as  indicated  by  the  shape  error  D  (cf.  Figure  2),  whereas  the  others 
apparently  did  not.  On  the  other  hand,  the  (CD-LF)  computations  appeared  to 
settle  into  a  permanent- form  solution  of  its  own,  a  property  not  evident  with  the 
other  schemes. 

In  purely  practical  terms  we  found  the  convergence  study  invaluable,  by  way 
of  exposing  errors  both  in  the  programming  and  of  a  conceptual  kind,  and  as  a 
guide  to  performing  the  computations  described  in  Bona,  Pritchard  and  Scott  (6], 
l7l,  and  (8). 

Finally,  the  numerical  evidence  that  the  (CD-LF)  scheme  might  have 
permanent- form  ( sol  it  ary- wave)  solution  was  a  surprise  to  us  and  it  would  be  of 
interest  to  know  definitively  whether  or  not  this  is  the  case.  (We  have  not  as 
yet  attempted  to  find  an  explicit  solitary-wave  solution  to  the  present  problem, 
or  to  demonstrate  the  existence  of  such  by  an  abstract  argument.)  Should  a 
family  of  such  solutions  exist  it  would  be  interesting  to  enquire  whether  or  not 
they  would  exhibit  the  so-called  soliton  property,  namely  that  two  solitary  waves 
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of  the  family  would  reemerge  from  an  interaction  with  their  shapes  unaltered.  We 
believe  the  answer  to  this  question  to  be  in  the  negative  for  the  following 
reason.  If  the  (CD-LF)  scheme  were  to  have  the  soliton  property  for  all  At  and 
Ax,  then  the  convergence  estimates  referred  to  in  §2  would  imply  that  the  same 
holds  for  (P):  as  At  and  Ax  tend  to  zero,  the  solitary-wave  solutions  to  (CD-LF) 
would  converge  to  a  solitary-wave  solution  of  (P).  But  the  numerical  experiments 
described  in  (6)  indicate  that  (P)  does  not  have  the  soliton  property.  It 
should,  however,  be  stressed  that  the  above  argument  is  not  a  proof,  even  given 
the  existence  of  ( CD-LF)  solitary  waves. 
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